
def LXM(s, m=2891336453, a=1310709051, x0=11565345, x1=5242836, k=64):
    
    # https://gist.github.com/trietptm/5cd60ed6add5adad6a34098ce255949a
    rotl = lambda val, r_bits, max_bits: \
        (val << r_bits%max_bits) & (2**max_bits-1) | \
        ((val & (2**max_bits-1)) >> (max_bits-(r_bits%max_bits)))    
    
    while True:
        # Combining operation
        z = s + x0
        # Mixing function (lea64)
        z = (z ^ (z >> 32)) * 0xdaba0b6eb09322e3
        z = (z ^ (z >> 32)) * 0xdaba0b6eb09322e3
        z = (z ^ (z >> 32))
        z &= (2**k) - 1
        # Update the LCG subgenerator
        s = (m * s + a)
        # Update the XBG subgenerator (xoroshiro128v1_0)
        q0 = x0; q1 = x1
        q1 ^= q0
        q0 = rotl(q0, 24, k)
        q0 = q0 ^ q1 ^ (q1 << 16)
        q1 = rotl(q1, 37, k)
        x0 = q0; x1 = q1
        # Return result (modified to be in interval 0-1)
        yield z/2**k


n = 5
print(f"{n} números generados al azar, en el intervalo 0-1:")

g = LXM(42) 
for r, i in zip(g,range(n)):
    print(f"* {r}")

5 números generados al azar, en el intervalo 0-1:
* 0.793781427453092
* 0.2155520740899848
* 0.7433470653358116
* 0.28660629734277526
* 0.515362051632894


scatterplot.show()


scatter3D.show()

Text(0.5, 0, 'X(n)+2')


histogramas.show()


boxplot.show()


qqplot.show()


from scipy.stats import chi2
import numpy as np
import math

class Chi2:
    def __init__(self, numero_de_muestras=1000, g = None):
        self._generador = g or LXM(42)
        self._numero_de_muestras = numero_de_muestras
        
    def realizar_prueba(self):
        frecuencias_observadas = self.frecuencias_observadas()
        
        f_esperada = 1/10 * sum(frecuencias_observadas)
        d2 = sum([(f_observada - f_esperada)**2 for f_observada in frecuencias_observadas]) / f_esperada

        limite_superior = chi2.ppf(0.95, df=5)
        print(f'Frecuencias: {frecuencias_observadas}')
        print('Estadístico: {:.2f}'.format(d2))
        print(f'Límite superior: {round(limite_superior, 2)}')
        self.validar_hipotesis(d2, limite_superior)

    def calcular_d2(self):
        numerador = (f_observada - f_esperada) ** 2
        denominador =  f_esperada
        d_al_cuadrado = numerador / denominados # acá es sumatoria
    
    def validar_hipotesis(self, d2, limite_superior):
        if d2 <= limite_superior:
            print("El test acepta la hipotesis nula")
        else:
            print("El test rechaza la hipótesis nula")
        
    def frecuencias_observadas(self):
        # Estas son las que tengo que obtener
        numeros_generados = []
        for i in range(self._numero_de_muestras):
            numeros_generados.append(next(self._generador))
        frecuencias = [0] * 10
        for numero in numeros_generados:
            pos = math.floor(numero*10)
            if pos < 10:
                frecuencias[pos] += 1
        return frecuencias


def pseudo_chi2():
    # Obtengo freuencias observadas de 0 a 9 (enteros)
    # La frecuencia esperada es 1/10  por la cantidad de números aleatorios que genero
    # Con estos valores ahora entonces calculo d2
    # d2 = sum (frecuencia_observada - frecuencia_esperada) ** 2  / frecuencia_esperada
    # Finalmente para chi2.ppf(0.95, df=5) vemos si d2 es mayor o menos al limite superior 
    return


# Corremos nuestro experimento
Chi2(numero_de_muestras=10000).realizar_prueba()

Frecuencias: [949, 998, 1019, 1007, 1008, 998, 1029, 968, 1001, 1023]
Estadístico: 5.48
Límite superior: 11.07
El test acepta la hipotesis nula


import matplotlib.pyplot as plt

class Gap:
    def __init__(self, cota_inferior, cota_superior, g = None):
        self._generador = g or LXM(42)
        self._cota_inferior = cota_inferior
        self._cota_superior = cota_superior
        self._p_0 = cota_superior - cota_inferior
        self._cantidad_de_gaps = 1000
        
    def realizar_prueba(self):
        # Espero a la primera vez que tengo valor dentro de la cota
        gaps = self.generar_gaps()
        self.plotear_gaps(gaps)
        frecuencias_observadas = self.generar_frecuencias(gaps)
        d2 = 0
        
        # Aplicamos chi2
        for i, frecuencia_observada in enumerate(frecuencias_observadas):
            f_esperada = self.frecuencia_esperada(i)
            numerador = (frecuencia_observada - f_esperada) ** 2
            denominador = f_esperada
            d2 += numerador / denominador
        
        # Validamos hipótesis
        k = len(frecuencias_observadas)
        limite_superior = chi2.ppf(0.95, df=k-1)
        print('Estadístico: {:.2f}'.format(d2))
        print(f'Límite superior: {round(limite_superior, 2)}')
        self.validar_hipotesis(d2, limite_superior)
        
    def validar_hipotesis(self, d2, limite_superior):
        if d2 <= limite_superior:
            print("El test acepta la hipotesis nula")
        else:
            print("El test rechaza la hipótesis nula")
            
    def generar_gaps(self):
        gaps = []
        # Generamos 1000 gaps
        while(len(gaps) < self._cantidad_de_gaps):
            while (True):
                n = next(self._generador)
                if n < self._cota_superior and n > self._cota_inferior:
                    break
            gap = 0
            while(True):
                n = next(self._generador)
                if n < self._cota_superior and n > self._cota_inferior:
                    gaps.append(gap)
                    break
                else:
                    gap +=1
        return gaps
    
    def plotear_gaps(self, gaps):
        plt.hist(gaps, color='b', bins="sturges", alpha=0.8, density=True, rwidth=0.8, align='left')
        plt.show()
        
    def generar_frecuencias(self, gaps):
        frecuencias = [0] * (max(gaps) + 1)
        for gap in gaps:
            frecuencias[gap] += 1
        return frecuencias
    
    def frecuencia_esperada(self, n):
        return (((1 - self._p_0)**n) * self._p_0) * self._cantidad_de_gaps


def pseudo_gap():
    # Obtengo gaps entre un rango de valores 
    # (dentro de [0,1]  ya que son los valores que genera el LXM programado
    # Con estos valores ahora entonces calculo d2
    # d2 = sum (frecuencia_observada - frecuencia_esperada) ** 2  / frecuencia_esperada
    # La  frecuencia esperada podrá ser distinta para cada valor, ya que
    # la misma sigue una distribución geométrica
    # Finalmente para chi2.ppf(0.95, df=cantidad_de_gaps) vemos si d2 es mayor o menos al limite superior 
    return


# Corremos nuestro experimento
Gap(cota_inferior=0.2, cota_superior=0.5).realizar_prueba()

Estadístico: 21.92
Límite superior: 31.41
El test acepta la hipotesis nula


from scipy.stats import chi2_contingency


class TestIndependencia:
    def __init__(self, g = None):
        self._generador = g or LXM(42)
        self._cantidad_de_vectores = 500
        
    def realizar_prueba(self):
        nros_x = []
        nros_y = []
        
        for i in range(self._cantidad_de_vectores):
            nros_x.append(next(self._generador))
            nros_y.append(next(self._generador))
        
        self.plot_vectores(nros_x, nros_y)
        
        d2, p, dof, ex = chi2_contingency([nros_x, nros_y], correction=True)
        
        # Validamos hipótesis
        limite_superior = chi2.ppf(0.95, df=dof)
        print('Estadístico: {:.2f}'.format(d2))
        print(f'Límite superior: {round(limite_superior, 2)}')
        self.validar_hipotesis(d2, limite_superior)
        
    def validar_hipotesis(self, d2, limite_superior):
        if d2 <= limite_superior:
            print("El test acepta la hipotesis nula")
        else:
            print("El test rechaza la hipótesis nula")

    def plot_vectores(self, nros_x, nros_y):
        plt.title(f"Dispersión de {self._cantidad_de_vectores} números generados al azar")
        plt.ylabel("X(n)-1")
        plt.xlabel("X(n)")
        plt.scatter(nros_x,nros_y,50)
        plt.show()


# Corremos nuestro experimento
TestIndependencia().realizar_prueba()

Estadístico: 87.60
Límite superior: 552.07
El test acepta la hipotesis nula


gen = LXM(42)

def generar_exponencial(lambda_):           
    u_1 = next(gen)
    return -np.log(u_1)/lambda_

def cociente(t):                            
    return e**((-(t-1)**2)/2)

def transformar_gaussiana(x, mean, std):    
    return x*std + mean

def generar_normal_0_1():
    no_aceptada = True
    while(no_aceptada):
        exp_1 = generar_exponencial(1)
        prob = cociente(exp_1)                
        if(next(gen) < prob):
            no_aceptada = False
            z = exp_1
            if(next(gen) < 0.5):
                z = -z
    return z


m1, m2, m3 = [], [], []     #muestras
for i in range(100):
    z1 = generar_normal_0_1()
    z2 = transformar_gaussiana(z1, 10, 2)
    m1.append(z2)

for i in range(1000):
    z1 = generar_normal_0_1()
    z2 = transformar_gaussiana(z1, 10, 2)
    m2.append(z2)

for i in range(10000):
    z1 = generar_normal_0_1()
    z2 = transformar_gaussiana(z1, 10, 2)
    m3.append(z2)


fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(20, 5))
fig.tight_layout(pad=3)
ax1.hist(m1, 50)
ax1.set_xlabel('Valores generados', fontsize=15)
ax1.set_title("Histograma (n = 100)", fontsize=20)
ax2.hist(m2, 50)
ax2.set_xlabel('Valores generados', fontsize=15)
ax2.set_title("Histograma (n = 1000)", fontsize=20)
ax3.hist(m3, 50)
ax3.set_xlabel('Valores generados', fontsize=15)
ax3.set_title("Histograma (n = 10000)", fontsize=20)

plt.show();


m1_ks, m2_ks, m3_ks = [], [], []     #muestras
for i in range(100):
    z1 = generar_normal_0_1()
    z2 = transformar_gaussiana(z1, 10, 2)
    m1_ks.append(z2)

for i in range(1000):
    z1 = generar_normal_0_1()
    z2 = transformar_gaussiana(z1, 10, 2)
    m2_ks.append(z2)

for i in range(5000):
    z1 = generar_normal_0_1()
    z2 = transformar_gaussiana(z1, 10, 2)
    m3_ks.append(z2)


state = np.random.default_rng()
dist_teorica_1 = sp.stats.norm(loc=10, scale=2).rvs(size=100, random_state=state)
dist_teorica_2 = sp.stats.norm(loc=10, scale=2).rvs(size=1000, random_state=state)
dist_teorica_3 = sp.stats.norm(loc=10, scale=2).rvs(size=5000, random_state=state)


figsts, (ax1sts, ax2sts, ax3sts) = plt.subplots(1, 3, figsize=(10, 3))
figsts.tight_layout(pad=3)

yval1 = np.linspace(0.0, 1.0,num=100)
ax1sts.step(x=np.sort(m1_ks), y=yval1)
ax1sts.step(x=np.sort(dist_teorica_1), y=yval1)
ax1sts.set_xlabel('Valores generados (n=100)', fontsize=8)
ax1sts.set_title("Probabilidad Acumulada", fontsize=10)

yval2 = np.linspace(0.0, 1.0,num=1000)
ax2sts.step(x=np.sort(m2_ks), y=yval2)
ax2sts.step(x=np.sort(dist_teorica_2), y=yval2)
ax2sts.set_xlabel('Valores generados (n=1000)', fontsize=8)
ax2sts.set_title("Probabilidad Acumulada", fontsize=10)

yval3 = np.linspace(0.0, 1.0,num=5000)
ax3sts.step(x=np.sort(m3_ks), y=yval3)
ax3sts.step(x=np.sort(dist_teorica_3), y=yval3)
ax3sts.set_xlabel('Valores generados (n=5000)', fontsize=8)
ax3sts.set_title("Probabilidad Acumulada", fontsize=10)

Text(0.5, 1.0, 'Probabilidad Acumulada')


figks, (ax1ks, ax2ks, ax3ks) = plt.subplots(1, 3, figsize=(20, 5))
figks.tight_layout(pad=3)
ax1ks.hist(m1_ks, 10)
ax1ks.set_xlabel('Valores generados', fontsize=15)
ax1ks.set_title("Histograma (n = 100)", fontsize=20)
ax2ks.hist(m2_ks, 25)
ax2ks.set_xlabel('Valores generados', fontsize=15)
ax2ks.set_title("Histograma (n = 1000)", fontsize=20)
ax3ks.hist(m3_ks, 50)
ax3ks.set_xlabel('Valores generados', fontsize=15)
ax3ks.set_title("Histograma (n = 5000)", fontsize=20)

plt.show();


from scipy import stats
ks1 = stats.kstest(m1_ks, dist_teorica_1, alternative='two-sided')
ks2 = stats.kstest(m2_ks, dist_teorica_2, alternative='two-sided')
ks3 = stats.kstest(m3_ks, dist_teorica_3, alternative='two-sided')

display ("Valor del estadístico del test Kolmogorov-Smirnov con n = 100: " + str(ks1.statistic))
display ("P-valor del test Kolmogorov-Smirnov con n = 100: " + str(ks1.pvalue))
display ("Valor del estadístico del test Kolmogorov-Smirnov con n = 1000: " + str(ks2.statistic))
display ("P-valor del test Kolmogorov-Smirnov con n = 1000: " + str(ks2.pvalue))
display ("Valor del estadístico del test Kolmogorov-Smirnov con n = 5000: " + str(ks3.statistic))
display ("P-valor del test Kolmogorov-Smirnov con n = 5000: " + str(ks3.pvalue))

'Valor del estadístico del test Kolmogorov-Smirnov con n = 100: 0.08'

'P-valor del test Kolmogorov-Smirnov con n = 100: 0.9084105017744525'

'Valor del estadístico del test Kolmogorov-Smirnov con n = 1000: 0.061'

'P-valor del test Kolmogorov-Smirnov con n = 1000: 0.04839715079181246'

'Valor del estadístico del test Kolmogorov-Smirnov con n = 5000: 0.0144'

'P-valor del test Kolmogorov-Smirnov con n = 5000: 0.6777877521935483'


m1_sw, m2_sw, m3_sw = [], [], []     #muestras
for i in range(100):
  z1 = generar_normal_0_1()
  z2 = transformar_gaussiana(z1, 10, 2)
  m1_sw.append(z2)

for i in range(1000):
  z1 = generar_normal_0_1()
  z2 = transformar_gaussiana(z1, 10, 2)
  m2_sw.append(z2)

for i in range(5000):
  z1 = generar_normal_0_1()
  z2 = transformar_gaussiana(z1, 10, 2)
  m3_sw.append(z2)


figsw, (ax1sw, ax2sw, ax3sw) = plt.subplots(1, 3, figsize=(20, 5))
figsw.tight_layout(pad=3)
ax1sw.hist(m1_sw, 10)
ax1sw.set_xlabel('Valores generados', fontsize=15)
ax1sw.set_title("Histograma (n = 100)", fontsize=20)
ax2sw.hist(m2_sw, 25)
ax2sw.set_xlabel('Valores generados', fontsize=15)
ax2sw.set_title("Histograma (n = 1000)", fontsize=20)
ax3sw.hist(m3_sw, 50)
ax3sw.set_xlabel('Valores generados', fontsize=15)
ax3sw.set_title("Histograma (n = 5000)", fontsize=20)

plt.show();


from scipy import stats
shapiro1 = stats.shapiro(m1_sw)
shapiro2 = stats.shapiro(m2_sw)
shapiro3 = stats.shapiro(m3_sw)

display ("Valor del estadístico del test Shapiro-Wilk con n = 100: " + str(shapiro1.statistic))
display ("P-valor del test Shapiro-Wilk con n = 100: " + str(shapiro1.pvalue))
display ("Valor del estadístico del test Shapiro-Wilk con n = 1000: " + str(shapiro2.statistic))
display ("P-valor del test Shapiro-Wilk con n = 1000: " + str(shapiro2.pvalue))
display ("Valor del estadístico del test Shapiro-Wilk con n = 5000: " + str(shapiro3.statistic))
display ("P-valor del test Shapiro-Wilk con n = 5000: " + str(shapiro3.pvalue))

'Valor del estadístico del test Shapiro-Wilk con n = 100: 0.9904636740684509'

'P-valor del test Shapiro-Wilk con n = 100: 0.7020654678344727'

'Valor del estadístico del test Shapiro-Wilk con n = 1000: 0.9989786744117737'

'P-valor del test Shapiro-Wilk con n = 1000: 0.8619867563247681'

'Valor del estadístico del test Shapiro-Wilk con n = 5000: 0.9995430111885071'

'P-valor del test Shapiro-Wilk con n = 5000: 0.28750666975975037'


import os

wd = os.getcwd()
if 'TP1' not in wd: wd += '/TP1'

arribos_real = []
with open(f"{wd}/tiempos_entre_arribos.txt") as f:
    for l in f.readlines(): arribos_real.append(float(l.strip()))

tasa = len(arribos_real)/sum(arribos_real)
print(f"Tengo {len(arribos_real)} arribos en {sum(arribos_real)} horas")
print(f"Mi tasa de arribos empírica es de {tasa} arribos por hora")

Tengo 10000 arribos en 1011.1994386169154 horas
Mi tasa de arribos empírica es de 9.889245996494681 arribos por hora


# Para generar eventos de Poisson:
#  1. Generar muestras de una distribución exponencial Z, 
#     con parametro Lambda: z1, z2, z3...
#  2. Definir tiempos de eventos conteo de Poisson: 
#     t0 = 0, t1 = t0 + z1, t2 = t1 + z2....

def iter_exponencial(lmbda, generador, n):
    """Genera una secuencia que siga una distribución exponencial, 
       dado un generador que sigue una distribución uniforme"""
    u = np.array([next(generador) for i in range(int(n))])
    x = -np.log(1-u)/lmbda
    return iter(x)

def generar_poisson(limite, lmbda, generador):
    expo = iter_exponencial(lmbda, generador, lmbda*limite*1.5)
    arribos = []
    t = 0
    while t <= limite:
        arribo = next(expo)
        arribos.append(arribo)
        t += arribo
    return arribos

limite = 24 * 31 # Un mes
arribos = generar_poisson(limite, tasa, LXM(42))

print(f"Simulé {len(arribos)} arribos en {sum(arribos)} horas")
print(f"Mi tasa de arribos simulada es de {len(arribos)/sum(arribos)} arribos por hora")

Simulé 7300 arribos en 744.064522714941 horas
Mi tasa de arribos simulada es de 9.810977109033201 arribos por hora


def comparar_poissones(poissones_con_label, n = 15):
    for arribos, label in poissones_con_label:
        plt.step(np.cumsum(arribos[:n]),range(n),where= 'post',label=label)
    plt.legend()
    plt.title(f"Comparación de los primeros {n} arribos")
    plt.show()
    
comparar_poissones([(arribos_real, "Poisson Empirica"), (arribos, "Poisson Simulada")])


limite = 96
poissones = []
for i in range(1000):
    poissones.append(generar_poisson(limite, tasa, LXM(i)))
    
comparar_poissones([(p, f"Poisson Generada #{i}") for i, p in enumerate(poissones[:5])])


import math

# PMF de poisson, a mano, para comparar con los valores teóricos
pmf = lambda x, lmbda: ((lmbda ** x) * (math.e ** -lmbda)) / (math.factorial(x))

print("Probabilidad que el primer vehículo arribe antes de los 10 minutos.")
diez_min = 1/6


# La probabilidad de que haya al menos un arribo en 10 minutos equivale a
#   1 - la probabilidad de que no haya ningun arribo en 10 minutos
print(f"* Probabilidad Teórica => {1 - pmf(0, tasa*diez_min)}")

favorable = []
for p in poissones:
    favorable.append(p[0] <= diez_min)

print(f"* Probabilidad Simulada => {sum(favorable)/len(favorable)}")

Probabilidad que el primer vehículo arribe antes de los 10 minutos.
* Probabilidad Teórica => 0.8076055651533732
* Probabilidad Simulada => 0.818


print("Probabilidad que el undécimo vehículo arribe después de los 60 minutos.")
sesenta_min = 1

# La probabilidad de que el vehiculo 11 llegué despues de 60 minutos equivale a
#   la sumatoria de que hayan llegado 0,1,2,...,10 vehiculos en 60 minutos
# Ojo con los off-by-one! `for i in range(11)` === de 0 a 10
print(f"* Probabilidad Teórica => {sum([pmf(i, tasa*sesenta_min) for i in range(11)])}")

favorable = []
for p in poissones:
    arribos_acumulados = np.cumsum(p)
    # El arribo en el indice [10] === El undecimo arribo
    favorable.append(arribos_acumulados[10] >= 1)

print(f"* Probabilidad Simulada => {sum(favorable)/len(favorable)}")

Probabilidad que el undécimo vehículo arribe después de los 60 minutos.
* Probabilidad Teórica => 0.596893339441658
* Probabilidad Simulada => 0.579


# Como nuestra pmf a mano da overflow al calcular un número grande como lambda a la 750,
#   importamos directamente el modulo de scipy
from scipy.stats import poisson

print("Probabilidad que arriben al menos 750 vehículos antes de las 72 horas.")
setentidos_horas = 72

# La probabilidad de que arriben al menos 750 vehiculos en 72 horas equivale a
#   1 - la sumatoria de que lleguen 0,1,...,749 en 72 horas
print(f"* Probabilidad Teórica => {1 - sum([poisson.pmf(i, tasa*setentidos_horas) for i in range(750)])}")

favorable = []
for p in poissones:
    arribos_acumulados = np.cumsum(p)
    favorable.append(arribos_acumulados[749] <= 72)

print(f"* Probabilidad Simulada => {sum(favorable)/len(favorable)}")

Probabilidad que arriben al menos 750 vehículos antes de las 72 horas.
* Probabilidad Teórica => 0.08097848535892749
* Probabilidad Simulada => 0.084


class EJ5:
    
    def __init__(self, gdf, lxm):
        self.gdf = gdf.copy()
        self.lxm = lxm
        self.puntos_en_barrios = {}
        for barrio in self.gdf['BARRIO'].tolist():
            self.puntos_en_barrios[barrio] = [[],[]]

    def emitir_barrios(self):
        display(self.gdf['BARRIO'])

    #metodo auxiliar para trasladar un barrio al primer cuadrante
    def trasladar_barrio(self, datos_barrio):
        extremos_barrio = datos_barrio.geometry.bounds

        signo1 = 1
        signo2 = 1

        if(extremos_barrio.iloc[0]['minx'] < 0):
            signo1 = -1
        if(extremos_barrio.iloc[0]['miny'] < 0):
            signo2 = -1

        return datos_barrio.translate(extremos_barrio.iloc[0]['minx']*signo1, extremos_barrio.iloc[0]['miny']*signo2)

    #Generar n puntos en un barrio
    def generar_en_barrio(self, barrio, cantidad):
        x, y = [], []

        #obtengo la geometria del barrio, la traslado al primer cuadrante
        #y busco el maximo entre sus extremos para formar un cuadrado que los contenga
        datos_barrio = self.gdf[self.gdf['BARRIO']==barrio]
        g_barrio = datos_barrio['geometry']
        g_barrio_t = self.trasladar_barrio(datos_barrio)
        extremos_b = datos_barrio.bounds.to_numpy()
        extremos_t = g_barrio_t.bounds.to_numpy()
        maximo_t = np.amax(extremos_t)

        #genero puntos dentro del cuadrado que obtuve, y los acepto si están dentro del barrio
        x = []
        y = []
        while len(x) < cantidad:
            z1 = next(self.lxm) * maximo_t
            z2 = next(self.lxm) * maximo_t
            ps = shapely.geometry.Point(z1, z2)
 
            if g_barrio_t.contains(ps).any():
                x.append(z1 + datos_barrio.geometry.bounds.minx)
                y.append(z2 + datos_barrio.geometry.bounds.miny)

        self.puntos_en_barrios[barrio] = [x, y]


    #Generar n puntos en todos los barrios
    def generar_en_todos(self, n):
        for barrio in self.gdf['BARRIO'].tolist():
            self.generar_en_barrio(barrio, n)

    #Graficar un barrio y los puntos generados en el mismo
    def graficar_en_barrio(self, barrio):
        x = self.puntos_en_barrios[barrio][0]
        y = self.puntos_en_barrios[barrio][1]
        title_string = "Puntos generados en el barrio de {} (n = {})".format(barrio, len(x))
        
        fig, ax = plt.subplots(figsize=(10, 10))
        fig.tight_layout(pad=3)
        ax.set_title(title_string, fontsize=20)

        datos_barrio = self.gdf.loc[self.gdf['BARRIO']==barrio]
        datos_barrio.plot(ax=ax, color='white', edgecolor='black')
        ax.scatter(x, y, s=3, color="red")
        ax.grid()

        return fig, ax

    #Graficar el mapa completo y los puntos generados en un barrio
    def graficar_mapa_completo_barrio(self, nombre_barrio):
        x = self.puntos_en_barrios[nombre_barrio][0]
        y = self.puntos_en_barrios[nombre_barrio][1]
        title_string = "Puntos generados en el barrio de {} (n = {})".format(nombre_barrio, len(x))
        
        fig, ax = plt.subplots(figsize=(10, 10))
        fig.tight_layout(pad=3)
        ax.set_title(title_string, fontsize=20)
        self.gdf.plot(ax=ax, color='white', edgecolor='black')
        ax.scatter(x, y, s=3, color="red")
        ax.grid()
        
        return fig, ax

    #Graficar el mapa completo y los puntos generados en todos los barrios (si hay)
    def graficar_mapa_completo_todos(self):
        title_string = "Puntos generados en todos los barrios."
        fig, ax = plt.subplots(figsize=(10, 10))
        fig.tight_layout(pad=3)
        ax.set_title(title_string, fontsize=20)

        self.gdf.plot(ax=ax, color='white', edgecolor='black')

        for b in self.gdf['BARRIO'].tolist():
            x = self.puntos_en_barrios[b][0]
            y = self.puntos_en_barrios[b][1]
            ax.scatter(x, y, s=3, color="red")
        ax.grid()

        return fig, ax

    def help(self):
        display("generar_en_barrio('barrio', cantidad)")
        display("graficar_mapa_completo_barrio('barrio')")
        display("graficar_mapa_completo_todos()")
        display("graficar_en_barrio('barrio')")
        display("generar_en_todos(cantidad)")


ej5test = EJ5(geodataframe, LXM(42))
ej5test.generar_en_barrio('PUERTO MADERO', 100)
ej5test.generar_en_todos(100)


ej5test.graficar_en_barrio('PUERTO MADERO')

(<Figure size 1000x1000 with 1 Axes>,
 <AxesSubplot: title={'center': 'Puntos generados en el barrio de PUERTO MADERO (n = 100)'}>)


ej5test.graficar_mapa_completo_barrio('PUERTO MADERO')

(<Figure size 1000x1000 with 1 Axes>,
 <AxesSubplot: title={'center': 'Puntos generados en el barrio de PUERTO MADERO (n = 100)'}>)


ej5test.graficar_mapa_completo_todos()

(<Figure size 1000x1000 with 1 Axes>,
 <AxesSubplot: title={'center': 'Puntos generados en todos los barrios.'}>)


# Ahora, rehagamos partes del trabajo, pero basados en un GCL en vez de en LXM
# Son los resultados distintos?
padrones = [97568, 100029, 100033, 104030]
semilla = int(np.mean(padrones))

def GCL(x, a = 1013904223, c = 1664525, m = 232):
    while True:
        x = ((a * x + c) % m) / m
        yield x


n = 5
g = GCL(semilla) 
print(f"{n} números generados al azar, en el intervalo 0-1:")
for r, i in zip(g,range(n)):
    print(f"* {r}")

5 números generados al azar, en el intervalo 0-1:
* 0.9224137931034483
* 0.29819411399035617
* 0.5019334719098848
* 0.895838293535956
* 0.008367325211393422


scatterplot.show()


Chi2(numero_de_muestras=10000, g=GCL(semilla)).realizar_prueba()

Frecuencias: [1005, 996, 1013, 1018, 1002, 992, 987, 1020, 994, 973]
Estadístico: 1.94
Límite superior: 11.07
El test acepta la hipotesis nula


limite = 24 * 31 # Un mes
arribos_gcl = generar_poisson(limite, tasa, GCL(semilla))

print(f"Simulé {len(arribos_gcl)} arribos en {sum(arribos_gcl)} horas")
print(f"Mi tasa de arribos simulada es de {len(arribos_gcl)/sum(arribos_gcl)} arribos por hora")
comparar_poissones([(arribos_real, "Poisson Empirica"), (arribos, "Poisson Simulada con LXM"), (arribos_gcl, "Poisson Simulada con GCL")])

Simulé 7394 arribos en 744.0969199691178 horas
Mi tasa de arribos simulada es de 9.936877578134409 arribos por hora


ej5test = EJ5(geodataframe, GCL(semilla))
ej5test.generar_en_barrio('PUERTO MADERO', 100)
ej5test.generar_en_todos(100)


ej5test.graficar_en_barrio('PUERTO MADERO')

(<Figure size 1000x1000 with 1 Axes>,
 <AxesSubplot: title={'center': 'Puntos generados en el barrio de PUERTO MADERO (n = 100)'}>)


ej5test.graficar_mapa_completo_barrio('PUERTO MADERO')

(<Figure size 1000x1000 with 1 Axes>,
 <AxesSubplot: title={'center': 'Puntos generados en el barrio de PUERTO MADERO (n = 100)'}>)


ej5test.graficar_mapa_completo_todos()

(<Figure size 1000x1000 with 1 Axes>,
 <AxesSubplot: title={'center': 'Puntos generados en todos los barrios.'}>)


import os
from bs4 import BeautifulSoup

wd = os.getcwd()
if 'TP2' not in wd: wd += '/TP2'

DIR = f'{wd}/paginas'
PAGINAS = os.listdir(DIR)

links = [] # Lista de tuplas de src->dst
palabras = {} # Dict de pagina: palabras

for f in PAGINAS:
    with open(os.path.join(DIR, f)) as file:
        soup = BeautifulSoup(file.read())
        palabras[f] = soup.get_text().split()

        for link in soup.find_all('a'):
            dst = link.get('href').replace('http://','')
            links.append((f,dst))

links

[('jonvoight.html', 'angelinajolie.html'),
 ('jonvoight.html', 'angelinajolie.html'),
 ('jonvoight.html', 'bradpitt.html'),
 ('bradpitt.html', 'jenniferaniston.html'),
 ('bradpitt.html', 'angelinajolie.html'),
 ('bradpitt.html', 'martinscorcese.html'),
 ('bradpitt.html', 'angelinajolie.html'),
 ('robertdeniro.html', 'martinscorcese.html'),
 ('angelinajolie.html', 'jonvoight.html'),
 ('angelinajolie.html', 'bradpitt.html')]


import networkx as nx

G = nx.DiGraph()
G.add_edges_from(links)
nx.draw_shell(G, with_labels=True, node_color='white', edge_color='grey')


matriz_ady = nx.to_pandas_adjacency(G)
matriz_ady


suma_filas = matriz_ady.sum(axis=1)
matriz_page_rank = matriz_ady.div(suma_filas, axis=0)
matriz_page_rank.fillna(1/len(PAGINAS), inplace=True)
matriz_page_rank


import numpy as np

# Hacemos N iteraciones de la matriz
iteraciones = np.linalg.matrix_power(matriz_page_rank, 50)

page_rank = {p: iteraciones[0][i] for i, p in enumerate(matriz_page_rank.index)}
page_rank

{'jonvoight.html': 0.16216216216216203,
 'angelinajolie.html': 0.21621621621621548,
 'bradpitt.html': 0.2432432432432427,
 'jenniferaniston.html': 0.13513513513513484,
 'martinscorcese.html': 0.18918918918918876,
 'robertdeniro.html': 0.05405405405405395}


import random

BUSQUEDAS = 20

todas_palabras = set()
for p in palabras.values():
    todas_palabras.update(p)

for p in random.choices(list(todas_palabras), k=BUSQUEDAS):
    matches = [f for f in PAGINAS if p in palabras[f]]
    print(f"'{p}' aparece en: {sorted(matches, key=lambda x: -page_rank[x])}")

'Twelve' aparece en: ['bradpitt.html']
'Jolie,' aparece en: ['bradpitt.html']
'(1999).' aparece en: ['angelinajolie.html']
'Derailed.' aparece en: ['jenniferaniston.html']
'(born' aparece en: ['bradpitt.html', 'angelinajolie.html', 'martinscorcese.html', 'jonvoight.html', 'jenniferaniston.html', 'robertdeniro.html']
'Break-Up,' aparece en: ['jenniferaniston.html']
'years,' aparece en: ['jonvoight.html']
'productions' aparece en: ['bradpitt.html']
'nomination.' aparece en: ['bradpitt.html', 'jonvoight.html']
'beautiful' aparece en: ['angelinajolie.html']
'named' aparece en: ['bradpitt.html']
'four' aparece en: ['bradpitt.html']
'Robert' aparece en: ['martinscorcese.html', 'robertdeniro.html']
'Part' aparece en: ['robertdeniro.html']
'Bull,' aparece en: ['robertdeniro.html']
'Mario' aparece en: ['robertdeniro.html']
'Leprechaun' aparece en: ['jenniferaniston.html']
'Taxi' aparece en: ['martinscorcese.html', 'robertdeniro.html']
'Came' aparece en: ['jenniferaniston.html']
'stock.' aparece en: ['martinscorcese.html']


class Metricas():
    def __init__(self, tasa_arribos, tasa_servicio):
        self.tasa_arribos = tasa_arribos
        self.tasa_servicio = tasa_servicio
        self.arribos = []
        self.salidas = []
        self.tiempo_en_cola = []
        self.tiempos_totales = []
        self.cantidad_en_cola = []
        self.tiempos_de_servicio = []
        self.tiempos_de_servicio_cola = []
        self.tasa_servicio_media_segun_largo_cola = []
        self.longitudMediaTeorica = 0
        self.longitudMediaEmpirica = 0
        self.tiempoMedioTeorico = 0
        self.tiempoMedioEmpirico = 0

    def calcular_metricas_paper(self):
        self.longitudMediaEmpirica = np.mean(self.cantidad_en_cola)
        self.longitudMediaTeorica =  (self.tasa_arribos/self.tasa_servicio)*np.e**(self.tasa_arribos/self.tasa_servicio)
        self.tiempoMedioEmpirico = np.mean(self.tiempos_totales)
        self.tiempoMedioTeorico = np.e**(self.tasa_arribos/self.tasa_servicio)/(self.tasa_servicio)
        #print(f'Tiempo de espera promedio obtenido de la simulación: {self.tiempoMedioEmpirico}')
        #print(f'Tiempo de espera promedio según el modelo: {self.tiempoMedioTeorico}')
        #print(f'Longitud de la cola promedio obtenida de la simulación: {self.longitudMediaEmpirica}')
        #print(f'Longitud de la cola promedio según el modelo: {self.longitudMediaTeorica}')


#generador para modelar la llegada de solicitudes
def solicitud(env, servicio, cant_arribos, tasa_arribos, tasa_servicio, metricas):
    solID = 0
    while solID < cant_arribos:
        t_solicitud_media = 1/tasa_arribos
        #t_solicitud = np.random.exponential(t_solicitud_media)
        t_solicitud = generar_exponencial(tasa_arribos)
        yield env.timeout(t_solicitud)
        metricas.arribos.append(env.now)
        solID += 1
        env.process(captura_imagen(env, servicio, solID, tasa_servicio, metricas))


#generador para modelar el servicio de captura de imágenes
def captura_imagen(env, servicio, solID, tasa_servicio, metricas):
    with servicio.request(priority=1) as req:
        m = len(servicio.queue)
        n = max(1, m)
        media_tiempos_de_servicio = 1/(tasa_servicio*n)
        tiempo_de_servicio = generar_exponencial(tasa_servicio)
        tiempo_de_servicio_cola = tiempo_de_servicio/n
        metricas.tasa_servicio_media_segun_largo_cola.append(1/tiempo_de_servicio_cola)
        metricas.tiempos_de_servicio.append(tiempo_de_servicio)
        metricas.tiempos_de_servicio_cola.append(tiempo_de_servicio_cola)
        metricas.cantidad_en_cola.append(len(servicio.queue))
        inicio_espera = env.now
        yield req
        fin_espera = env.now
        inicio_servicio = env.now
        yield env.timeout(tiempo_de_servicio_cola)
        fin_servicio = env.now
        metricas.tiempo_en_cola.append(fin_espera - inicio_espera)
        metricas.tiempos_totales.append(fin_servicio - inicio_espera)
        metricas.salidas.append(env.now)


#inicio la simulación
def simular(cant_satelites, cant_arribos, tasa_arribos, tasa_servicio, semilla):

    metricas = Metricas(tasa_arribos, tasa_servicio)
    env = simpy.Environment()
    servicio = simpy.PriorityResource(env, capacity = cant_satelites)
    np.random.seed(seed=semilla)
    print(f'Simulando {cant_arribos} arribos con lambdaArribos = {tasa_arribos} y lambdaServicio = {tasa_servicio}')
    env.process(solicitud(env, servicio, cant_arribos, tasa_arribos, tasa_servicio, metricas))
    env.run()
    return metricas


metricasHist = []
a_arribos = np.linspace(0.5, 1.5, 3)
a_llegadas = np.linspace(0.5, 1.5, 3)
longitudes = []
tiempos = []
cantidades_y_tasa_media = []
k=0
for i in a_arribos:
    for j in a_llegadas:
        print(k)
        k=k+1
        metricas = simular(cant_satelites=1, cant_arribos=100, tasa_arribos=i, tasa_servicio=j, semilla=97568)
        metricas.calcular_metricas_paper()
        longitudes.append((metricas.longitudMediaTeorica, metricas.longitudMediaEmpirica))
        tiempos.append((metricas.tiempoMedioTeorico, metricas.tiempoMedioEmpirico))
        cantidades_y_tasa_media.append((np.mean(metricas.cantidad_en_cola), (np.mean(metricas.tasa_servicio_media_segun_largo_cola))))
        metricasHist.append(metricas)

0
Simulando 100 arribos con lambdaArribos = 0.5 y lambdaServicio = 0.5
1
Simulando 100 arribos con lambdaArribos = 0.5 y lambdaServicio = 1.0
2
Simulando 100 arribos con lambdaArribos = 0.5 y lambdaServicio = 1.5
3
Simulando 100 arribos con lambdaArribos = 1.0 y lambdaServicio = 0.5
4
Simulando 100 arribos con lambdaArribos = 1.0 y lambdaServicio = 1.0
5
Simulando 100 arribos con lambdaArribos = 1.0 y lambdaServicio = 1.5
6
Simulando 100 arribos con lambdaArribos = 1.5 y lambdaServicio = 0.5
7
Simulando 100 arribos con lambdaArribos = 1.5 y lambdaServicio = 1.0
8
Simulando 100 arribos con lambdaArribos = 1.5 y lambdaServicio = 1.5


fig, ax = plt.subplots(nrows=1, ncols=3, figsize=(20, 5))
ax[0].hist(metricasHist[3].cantidad_en_cola, 20,rwidth=0.9)
ax[1].hist(metricasHist[4].cantidad_en_cola, 20,rwidth=0.9)
ax[2].hist(metricasHist[5].cantidad_en_cola, 20,rwidth=0.9)
ax[0].set_xticks(np.arange(0, 20, 1))
ax[1].set_xticks(np.arange(0, 20, 1))
ax[2].set_xticks(np.arange(0, 20, 1))
ax[0].set_xlabel(f'Cantidad de solicitudes en cola \n (λ = {metricasHist[3].tasa_arribos}, μ = {metricasHist[3].tasa_servicio})', fontsize=15)
ax[1].set_xlabel(f'Cantidad de solicitudes en cola \n (λ = {metricasHist[4].tasa_arribos}, μ = {metricasHist[4].tasa_servicio})', fontsize=15)
ax[2].set_xlabel(f'Cantidad de solicitudes en cola \n (λ = {metricasHist[5].tasa_arribos}, μ = {metricasHist[5].tasa_servicio})', fontsize=15)
fig.suptitle("Histogramas de solicitudes en cola", fontsize=20)

Text(0.5, 0.98, 'Histogramas de solicitudes en cola')


fig, ax3 = plt.subplots(nrows=1, ncols=3, figsize=(20, 5))
ax3[0].step(x=metricasHist[3].arribos, y=range(len(metricasHist[3].arribos)))
ax3[0].step(x=metricasHist[3].salidas, y=range(len(metricasHist[3].salidas)))
ax3[1].step(x=metricasHist[4].arribos, y=range(len(metricasHist[4].arribos)))
ax3[1].step(x=metricasHist[4].salidas, y=range(len(metricasHist[4].salidas)))
ax3[2].step(x=metricasHist[5].arribos, y=range(len(metricasHist[5].arribos)))
ax3[2].step(x=metricasHist[5].salidas, y=range(len(metricasHist[5].salidas)))
fig.suptitle("Evolución de arribos y salidas", fontsize=20)
ax3[0].set_xlabel(f'Arribos y salidas \n (λ = {metricasHist[3].tasa_arribos}, μ = {metricasHist[3].tasa_servicio})', fontsize=15)
ax3[1].set_xlabel(f'Arribos y salidas \n (λ = {metricasHist[4].tasa_arribos}, μ = {metricasHist[4].tasa_servicio})', fontsize=15)
ax3[2].set_xlabel(f'Arribos y salidas \n (λ = {metricasHist[5].tasa_arribos}, μ = {metricasHist[5].tasa_servicio})', fontsize=15)
ax3[0].legend([f'λ = {metricasHist[3].tasa_arribos}, μ = {metricasHist[3].tasa_servicio}'])
ax3[1].legend([f'λ = {metricasHist[4].tasa_arribos}, μ = {metricasHist[4].tasa_servicio}'])
ax3[2].legend([f'λ = {metricasHist[5].tasa_arribos}, μ = {metricasHist[5].tasa_servicio}'])

<matplotlib.legend.Legend at 0x7f77a89bb8e0>


fig, ax4 = plt.subplots(nrows=2, ncols=1, figsize=(20, 10))
p = sns.kdeplot(metricasHist[0].cantidad_en_cola,ax=ax4[0])
p = sns.kdeplot(metricasHist[3].cantidad_en_cola,ax=ax4[0])
p = sns.kdeplot(metricasHist[6].cantidad_en_cola,ax=ax4[0])
#p = sns.kdeplot(metricasHist[12].cantidad_en_cola,ax=ax4[0])
p = sns.kdeplot(metricasHist[0].cantidad_en_cola,ax=ax4[1])
p = sns.kdeplot(metricasHist[1].cantidad_en_cola,ax=ax4[1])
p = sns.kdeplot(metricasHist[2].cantidad_en_cola,ax=ax4[1])
#p = sns.kdeplot(metricasHist[3].cantidad_en_cola,ax=ax4[1])
ax4[0].set_xticks(np.arange(0, 30, 1))
ax4[0].legend([f'λ = {metricasHist[0].tasa_arribos}, μ = {metricasHist[0].tasa_servicio}',
               f'λ = {metricasHist[3].tasa_arribos}, μ = {metricasHist[3].tasa_servicio}',
               f'λ = {metricasHist[6].tasa_arribos}, μ = {metricasHist[6].tasa_servicio}'])
ax4[0].set_xlabel('Solicitudes en cola', fontsize=15)
ax4[0].set_ylabel('Densidad', fontsize=15)
ax4[0].set_title("Densidad de cantidad de solicitudes en cola", fontsize=20)
ax4[1].set_xticks(np.arange(0, 30, 1))
ax4[1].legend([f'λ = {metricasHist[0].tasa_arribos}, μ = {metricasHist[0].tasa_servicio}',
               f'λ = {metricasHist[1].tasa_arribos}, μ = {metricasHist[1].tasa_servicio}',
               f'λ = {metricasHist[2].tasa_arribos}, μ = {metricasHist[2].tasa_servicio}'])
ax4[1].set_xlabel('Solicitudes en cola', fontsize=15)
ax4[1].set_ylabel('Densidad', fontsize=15)
ax4[1].set_title("Densidad de cantidad de solicitudes en cola", fontsize=20)
fig.tight_layout()


fig, ax5 = plt.subplots(nrows=1, ncols=1, figsize=(20, 5))
ax5.scatter(*zip(*cantidades_y_tasa_media))
ax5.set_xlabel('Mean image capture queue length', fontsize=15)
ax5.set_ylabel('Mean image capture service rate (req/day)', fontsize=15)
ax5.set_title("Longitudes de cola medias", fontsize=20)

Text(0.5, 1.0, 'Longitudes de cola medias')


fig, ax6 = plt.subplots(nrows=1, ncols=2, figsize=(20, 5))
ax6[0].scatter(*zip(*longitudes))
ax6[0].set_xlabel('Longitud Media Teórica', fontsize=15)
ax6[0].set_ylabel('Longitud Media Empirica', fontsize=15)
ax6[0].set_title("Longitudes de cola medias", fontsize=20)
ax6[1].scatter(*zip(*tiempos))
ax6[1].set_xlabel('Tiempo medio de servicio Teórico', fontsize=15)
ax6[1].set_ylabel('Tiempo medio de servicio Empírico', fontsize=15)
ax6[1].set_title("Tiempos medios de servicio", fontsize=20)

Text(0.5, 1.0, 'Tiempos medios de servicio')


fig, ax7 = plt.subplots(nrows=3, ncols=2, figsize=(20,10))
p = sns.kdeplot(metricasHist[3].tiempos_de_servicio,ax=ax7[0,0])
p = sns.kdeplot(metricasHist[3].tiempos_de_servicio_cola,ax=ax7[0,1])
p = sns.kdeplot(metricasHist[4].tiempos_de_servicio,ax=ax7[1,0])
p = sns.kdeplot(metricasHist[4].tiempos_de_servicio_cola,ax=ax7[1,1])
p = sns.kdeplot(metricasHist[5].tiempos_de_servicio,ax=ax7[2,0])
p = sns.kdeplot(metricasHist[5].tiempos_de_servicio_cola,ax=ax7[2,1])
ax7[0,0].legend([f'λ = {metricasHist[3].tasa_arribos}, μ = {metricasHist[3].tasa_servicio}'])
ax7[0,1].legend([f'λ = {metricasHist[3].tasa_arribos}, μ = {metricasHist[3].tasa_servicio}'])
ax7[1,0].legend([f'λ = {metricasHist[4].tasa_arribos}, μ = {metricasHist[4].tasa_servicio}'])
ax7[1,1].legend([f'λ = {metricasHist[4].tasa_arribos}, μ = {metricasHist[4].tasa_servicio}'])
ax7[2,0].legend([f'λ = {metricasHist[5].tasa_arribos}, μ = {metricasHist[5].tasa_servicio}'])
ax7[2,1].legend([f'λ = {metricasHist[5].tasa_arribos}, μ = {metricasHist[5].tasa_servicio}'])
ax7[0,0].set_yticks(np.arange(0, 1.6, 0.2))
ax7[0,1].set_yticks(np.arange(0, 1.6, 0.2))
ax7[1,0].set_yticks(np.arange(0, 1.6, 0.2))
ax7[1,1].set_yticks(np.arange(0, 1.6, 0.2))
ax7[2,0].set_yticks(np.arange(0, 1.6, 0.2))
ax7[2,1].set_yticks(np.arange(0, 1.6, 0.2))
ax7[0,1].legend([f'λ = {metricasHist[3].tasa_arribos}, μ = {metricasHist[3].tasa_servicio}'])
ax7[1,0].legend([f'λ = {metricasHist[4].tasa_arribos}, μ = {metricasHist[4].tasa_servicio}'])
ax7[1,1].legend([f'λ = {metricasHist[4].tasa_arribos}, μ = {metricasHist[4].tasa_servicio}'])
ax7[2,0].legend([f'λ = {metricasHist[5].tasa_arribos}, μ = {metricasHist[5].tasa_servicio}'])
ax7[2,1].legend([f'λ = {metricasHist[5].tasa_arribos}, μ = {metricasHist[5].tasa_servicio}'])
#fig.tight_layout()
fig.suptitle("Densidad de tiempos de servicio", fontsize=20)

Text(0.5, 0.98, 'Densidad de tiempos de servicio')


from scipy.stats import expon, uniform
import seaborn as sns
import matplotlib.pyplot as plt
import matplotlib.pylab as pylab
from math import floor

sns.set(style="darkgrid")
params = {
    'axes.titlesize': 20,
    'axes.labelsize': 15,
    'legend.fontsize': 20,
    'figure.figsize': (10, 6),
}
pylab.rcParams.update(params)


class Consulta:
    def __init__(self, tiempo_arribo, duracion, tiempo_finalizacion_anterior):
        if tiempo_arribo < tiempo_finalizacion_anterior:
            self.duracion_espera = tiempo_finalizacion_anterior - tiempo_arribo
            self.tuvo_que_esperar = True
        else:
            self.duracion_espera = 0
            self.tuvo_que_esperar = False
                            
        self.tiempo_finalizacion = tiempo_arribo + duracion + self.duracion_espera


def simular_consultas(tiempos_de_arribo, duraciones):
    consultas = []

    # Uso -1 como tiempo de finalizacion del anterior para que empieze ni bien llega
    primera_consulta = Consulta(tiempos_de_arribo[0], duraciones[0], tiempo_finalizacion_anterior = -1)
    consultas.append(primera_consulta)

    for tiempo_arribo, duracion in zip(tiempos_de_arribo[1:], duraciones[1:]):
        consulta = Consulta(tiempo_arribo,
                        duracion, 
                        tiempo_finalizacion_anterior = consultas[-1].tiempo_finalizacion)
        consultas.append(consulta)
        
    tiempo_espera_total = 0
    cant_consultas_que_no_esperaron = 0
    for consulta in consultas:
        tiempo_espera_total += consulta.duracion_espera
        if not consulta.tuvo_que_esperar:
            cant_consultas_que_no_esperaron += 1

    return tiempo_espera_total, cant_consultas_que_no_esperaron


def mostrar_resultados(titulo, tiempo_espera_total, cant_consultas_que_no_esperaron, cant_consultas):
    print(f"--------------{titulo}--------------")
    
    
    segundos = round(tiempo_espera_total / cant_consultas, 0)
    
    minutos = floor(segundos/60)
    if(minutos != 0):
        segundos = segundos % (minutos * 60)
    
    horas = floor(minutos/60)
    if(horas != 0):
        minutos = minutos % (horas * 60)

        
    dias = floor(horas/24)
    if(dias != 0):
        horas = horas % (dias * 24)
    
    print(f"\tEl tiempo de espera promedio fue ", end="")
    if(dias != 0):
        print(str(dias) + " Dias ", end="")
    if(horas != 0):    
        print(str(horas) + " Horas ", end="")
    if(minutos != 0):    
        print(str(minutos) + " Minutos ", end="")
    if(segundos != 0):    
        print(str(segundos) + " Segundos")
    
    
    print(f"\tEl {round(cant_consultas_que_no_esperaron/cant_consultas, 5)}% de las consultas no tuvo que esperar")


def simular_consultas_centralizado(cant_consultas, media_tiempos_de_arribo, media_duraciones_consultas):
    tiempos_de_arribo = expon.rvs(size=cant_consultas, scale=media_tiempos_de_arribo, random_state=10)
    duraciones = expon.rvs(size=cant_consultas, scale=media_duraciones_consultas, random_state=10)
    return simular_consultas(tiempos_de_arribo, duraciones)


t, f = simular_consultas_centralizado(100000, 4, 0.8)
mostrar_resultados("Base de datos central", t,f, 100000)

--------------Base de datos central--------------
	El tiempo de espera promedio fue 11 Horas 3 Minutos 16.0 Segundos
	El 2e-05% de las consultas no tuvo que esperar


def simular_consultas_distribuido(cant_consultas, media_tiempos_de_arribo,
                                  media_duraciones_consultas1, media_duraciones_consultas2, p):
    tiempos_de_arribo = expon.rvs(size=cant_consultas, scale=media_tiempos_de_arribo, random_state=10)
    u = uniform.rvs(size=cant_consultas, random_state=10)
    
    tiempos_de_arribo1 = []
    tiempos_de_arribo2 = []
    
    for ui, tiempo_de_arribo in zip(u, tiempos_de_arribo):
        if p <= ui:
            tiempos_de_arribo1.append(tiempo_de_arribo)
        else:
            tiempos_de_arribo2.append(tiempo_de_arribo)
    
    
    duraciones1 = expon.rvs(size=len(tiempos_de_arribo1), scale=media_duraciones_consultas1, random_state=10)
    t1, f1 = simular_consultas(tiempos_de_arribo1, duraciones1)
    
    duraciones2 = expon.rvs(size=len(tiempos_de_arribo2), scale=media_duraciones_consultas2, random_state=10)
    t2, f2 = simular_consultas(tiempos_de_arribo2, duraciones2)
    
    return t1, f1, len(tiempos_de_arribo1), t2, f2, len(tiempos_de_arribo2)


t1, f1, l1, t2, f2, l2 = simular_consultas_distribuido(100000, 4, 0.7, 1, 0.7)
mostrar_resultados("Ambas bases distribuidas", t1 + t2, f2 + f2, l1 + l2)
mostrar_resultados("Base de datos distribuida 1", t1, f1, l1)
mostrar_resultados("Base de datos distribuida 2", t2, f2, l2)

--------------Ambas bases distribuidas--------------
	El tiempo de espera promedio fue 7 Horas 40 Minutos 21.0 Segundos
	El 4e-05% de las consultas no tuvo que esperar
--------------Base de datos distribuida 1--------------
	El tiempo de espera promedio fue 2 Horas 51 Minutos 34.0 Segundos
	El 7e-05% de las consultas no tuvo que esperar
--------------Base de datos distribuida 2--------------
	El tiempo de espera promedio fue 9 Horas 42 Minutos 30.0 Segundos
	El 3e-05% de las consultas no tuvo que esperar


def simular_variaciones(cantidades_consultas=[100000], 
                        medias_tiempos_de_arribo=[4],
                        medias_duraciones_consultas = [0.8],  
                        medias_duraciones_consultas1 = [0.7],
                        medias_duraciones_consultas2 = [1],
                        lista_p = [0.7]):

    tiempos_centralizado = []
    frecuencias_centralizado = []

    tiempos_distribuido = []
    frecuencias_distribuido = []
                         
    for cant_consultas in cantidades_consultas:
        for media_tiempos_de_arribo in medias_tiempos_de_arribo:
            for media_duraciones_consultas in medias_duraciones_consultas:
                for media_duraciones_consultas1 in medias_duraciones_consultas1:
                     for media_duraciones_consultas2 in medias_duraciones_consultas2:
                            for p in lista_p:
                                t, f = simular_consultas_centralizado(cant_consultas, media_tiempos_de_arribo, 
                                                                      media_duraciones_consultas)
                                tiempos_centralizado.append(t / cant_consultas)
                                frecuencias_centralizado.append(f /cant_consultas)
    
                                t1, f1, l1, t2, f2, l2  = simular_consultas_distribuido(cant_consultas, 
                                                                                        media_tiempos_de_arribo,
                                                                                        media_duraciones_consultas1,
                                                                                        media_duraciones_consultas2,
                                                                                        p)
    
                                tiempos_distribuido.append((t1 + t2) / (l1 + l2))
                                frecuencias_distribuido.append((f1 + f2) / (l1 + l2))
            
            
    return tiempos_centralizado, frecuencias_centralizado, tiempos_distribuido, frecuencias_distribuido


def graficar_tiempo_espera(x, xlabel, tiempos_centralizado, tiempos_distribuido,
                         centralizado_es_referencia = False, distribuido_es_referencia = False):
    
    if centralizado_es_referencia:
        color_centralizado = 'b'
        linestyle_centralizado = '--'
    else:
        color_centralizado = 'bo'
        linestyle_centralizado = '-'
        
    if distribuido_es_referencia:
        color_distribuido = 'g'
        linestyle_distribuido = '--'
    else:
        color_distribuido = 'go'
        linestyle_distribuido = '-'
    
    plt.plot(x, tiempos_centralizado, color_centralizado, linestyle=linestyle_centralizado, label="Base centralizada")
    plt.plot(x, tiempos_distribuido, color_distribuido, linestyle=linestyle_distribuido, label="Bases distribuidas")
    plt.title(f"Tiempo de espera medio según\n{xlabel.lower()}")
    plt.xlabel(xlabel)
    plt.ylabel("Tiempo de espera medio")
    plt.legend()
    plt.show()


def graficar_fraccion_que_no_espero(x, xlabel, tiempos_centralizado, tiempos_distribuido,
    centralizado_es_referencia = False, distribuido_es_referencia = False):
    
    if centralizado_es_referencia:
        color_centralizado = 'b'
        linestyle_centralizado = '--'
    else:
        color_centralizado = 'bo'
        linestyle_centralizado = '-'
        
    if distribuido_es_referencia:
        color_distribuido = 'g'
        linestyle_distribuido = '--'
    else:
        color_distribuido = 'go'
        linestyle_distribuido = '-'
    
    plt.plot(x, tiempos_centralizado, color_centralizado, linestyle=linestyle_centralizado, label="Base centralizada")
    plt.plot(x, tiempos_distribuido, color_distribuido, linestyle=linestyle_distribuido, label="Bases distribuidas")
    plt.title(f"Fracción de solicitudes que no esperaron según\n{xlabel.lower()}")
    plt.xlabel(xlabel)
    plt.ylabel("Fracción de solicitudes que no esperaron")
    plt.legend()


medias_tiempos_de_arribo = range(1, 11)
tiempos_centralizado, frecuencias_centralizado, tiempos_distribuido, frecuencias_distribuido = simular_variaciones(
    medias_tiempos_de_arribo = medias_tiempos_de_arribo
)


graficar_tiempo_espera(medias_tiempos_de_arribo, "Media de los tiempos de arribo de consultas",
                      tiempos_centralizado, tiempos_distribuido)


graficar_fraccion_que_no_espero(medias_tiempos_de_arribo, "Media de los tiempos de arribo de consultas",
                      frecuencias_centralizado, frecuencias_distribuido)


lambda_l = []
for i in range(2, 20, 2):
    lambda_l.append(i/10)
tiempos_centralizado, frecuencias_centralizado, tiempos_distribuido, frecuencias_distribuido = simular_variaciones(
     medias_duraciones_consultas = lambda_l
)


graficar_tiempo_espera(lambda_l, "Lambda",
                      tiempos_centralizado, tiempos_distribuido, distribuido_es_referencia = True)


graficar_fraccion_que_no_espero(lambda_l, "Lambda",
                      frecuencias_centralizado, frecuencias_distribuido, distribuido_es_referencia = True)


lambda_1 = []
for i in range(2, 20, 2):
    lambda_1.append(i/10)
tiempos_centralizado, frecuencias_centralizado, tiempos_distribuido, frecuencias_distribuido = simular_variaciones(
     medias_duraciones_consultas1 = lambda_1
)


graficar_tiempo_espera(lambda_1, "Lambda1",
                      tiempos_centralizado, tiempos_distribuido, centralizado_es_referencia = True)


graficar_fraccion_que_no_espero(lambda_1, "Lambda1",
                      frecuencias_centralizado, frecuencias_distribuido, centralizado_es_referencia = True)


lambda_2 = []
for i in range(2, 20, 2):
    lambda_2.append(i/10)
tiempos_centralizado, frecuencias_centralizado, tiempos_distribuido, frecuencias_distribuido = simular_variaciones(
     medias_duraciones_consultas1 = lambda_2
)


graficar_tiempo_espera(lambda_2, "Lambda2",
                      tiempos_centralizado, tiempos_distribuido, centralizado_es_referencia = True)


graficar_fraccion_que_no_espero(lambda_2, "Lambda2",
                      frecuencias_centralizado, frecuencias_distribuido, centralizado_es_referencia = True)


lista_p = []
for i in range(1, 10):
    lista_p.append(i/10)
tiempos_centralizado, frecuencias_centralizado, tiempos_distribuido, frecuencias_distribuido = simular_variaciones(
     lista_p = lista_p
)


graficar_tiempo_espera(lista_p, "p",
                      tiempos_centralizado, tiempos_distribuido, centralizado_es_referencia = True)


graficar_fraccion_que_no_espero(lista_p, "p",
                      frecuencias_centralizado, frecuencias_distribuido, centralizado_es_referencia = True)


from scipy.stats import expon, uniform
import matplotlib.pyplot as plt
import matplotlib.pylab as pylab
from math import floor
from random import randint
import random
import matplotlib.pyplot as plt
from statistics import median


class Cajero:
    def __init__(self):
        self._billetes_en_cajero = 2000 # son todos de 100,  arranca el  día lleno
        self._capacidad_maxima = 2000  #billetes max q entran
        self._historial_dia_actual = []
        self._historial = []
        
    def nuevo_dia(self):
        self._billetes_en_cajero = 2000
        self._historial.append(self._historial_dia_actual)
        self._historial_dia_actual = []
    
    def retirar(self, cantidad_de_billetes):
        if self._billetes_en_cajero < cantidad_de_billetes:
            return False
        self._billetes_en_cajero -= cantidad_de_billetes
        self._historial_dia_actual.append(self._billetes_en_cajero)
        return True
        
    def depositar(self, cantidad_de_billetes):
        if self._billetes_en_cajero + cantidad_de_billetes > self._capacidad_maxima:
            return False
        self._billetes_en_cajero += cantidad_de_billetes
        self._historial_dia_actual.append(self._billetes_en_cajero)
        return True
        
    def historial(self):
        return self._historial
        
        
class Banco:
    def __init__(self):
        self._cajero = Cajero()
        self._minutos = 0 # al pasa 24 * 60 termina el día y restarteamos
        self._tiempos_de_espera = []
        
    def _simular_dia(self):
        print(expon.rvs(size=3, random_state=10))
        # tarda en llegar n 
        
    def _grupo_cliente(self):
        return random.choices([1,2], weights=(75, 25))[0]
    
    def historial_de_cajero(self):
        return self._cajero.historial()
            
    def simular(self, dias, imprimir_estadisticas_clientes=False, imprimir_media_espera=False):
        tiempos_de_arribo = expon.rvs(size=144 * 60 * dias, scale=10, random_state=10)
        duraciones_1 = expon.rvs(size=len(tiempos_de_arribo), scale=1.5, random_state=10)
        duraciones_2 = expon.rvs(size=len(tiempos_de_arribo), scale=5, random_state=10)
        clientes_perdidos = 0
        cantidad_de_clientes = 0
        cantidad_de_dias = 0
        cantidad_de_clientes = 0
        minutos_en_hora = 24 * 60
        for i, tiempo in enumerate(tiempos_de_arribo):
            grupo_cliente =  self._grupo_cliente()
            self._minutos += tiempo
            duracion_de_persona = 0
            if self._minutos >= minutos_en_hora:
                self._minutos = minutos_en_hora - self._minutos
                self._cajero.nuevo_dia()
                cantidad_de_dias += 1
                if cantidad_de_dias == dias:
                    break
            if grupo_cliente == 1: # Retiran
                cant_billetes = randint(3, 50)
                fue_exitoso = self._cajero.retirar(cant_billetes)
                clientes_perdidos += 0 if fue_exitoso else 1
                duracion_de_persona = duraciones_1[i]
            else: # Depositan
                cant_billetes = randint(10, 110)
                fue_exitoso = self._cajero.depositar(cant_billetes)
                clientes_perdidos += 0 if fue_exitoso else 1
                duracion_de_persona = duraciones_2[i]
            demora = 0
            if i - 1 < len(tiempos_de_arribo)  and tiempos_de_arribo[i+1] < duracion_de_persona:
                demora = duracion_de_persona - tiempos_de_arribo[i+1]
                self._minutos += demora
            self._tiempos_de_espera.append(duracion_de_persona + demora)
            cantidad_de_clientes += 1
        if imprimir_estadisticas_clientes:
            print(f'Pasaron {cantidad_de_dias} días')
            print(f'Cantidad de clientes: {cantidad_de_clientes}')
            print(f'Se perdieron {clientes_perdidos} clientes')
            porcentaje_clientes_perdidos = clientes_perdidos * 100 / cantidad_de_clientes
            print(f'Porcentaje de clientes perdidos: {floor(porcentaje_clientes_perdidos)}%')
        if imprimir_media_espera:
            mediana = round(median(self._tiempos_de_espera), 2)
            print(f'Tiempo medio de espera: {mediana} minutos')


# Simular 1000 días
banco = Banco()
banco.simular(dias=1000, imprimir_estadisticas_clientes=True)

Pasaron 1000 días
Cantidad de clientes: 138501
Se perdieron 2587 clientes
Porcentaje de clientes perdidos: 1%


# cantidad de billetes en el cajero dps de cada transacción
banco = Banco()
banco.simular(dias=1000)
monto_en_cajero = banco.historial_de_cajero()[1]

# make up some data
x = [i for i,monto in enumerate(monto_en_cajero)]

# plot
plt.scatter(x,monto_en_cajero, linewidths = 0.1)
# beautify the x-labels
plt.gcf().autofmt_xdate()


#  Calcular el tiempo medio que los clientes demoraron en el sistema (espera + utilización del cajero)
banco = Banco()
banco.simular(dias=1000, imprimir_media_espera=True)

Tiempo medio de espera: 1.38 minutos

	jonvoight.html	angelinajolie.html	bradpitt.html	jenniferaniston.html	martinscorcese.html	robertdeniro.html
jonvoight.html	0.000000	0.500000	0.500000	0.000000	0.000000	0.000000
angelinajolie.html	0.500000	0.000000	0.500000	0.000000	0.000000	0.000000
bradpitt.html	0.000000	0.333333	0.000000	0.333333	0.333333	0.000000
jenniferaniston.html	0.166667	0.166667	0.166667	0.166667	0.166667	0.166667
martinscorcese.html	0.166667	0.166667	0.166667	0.166667	0.166667	0.166667
robertdeniro.html	0.000000	0.000000	0.000000	0.000000	1.000000	0.000000

Trabajos Prácticos de Simulación [75.26]¶

Trabajo Práctico 1¶

LXM¶

Test Estadísticos¶

Test de Chi2¶

Test de Gap¶

Test de Independencia¶

Ejercicio 3¶

Test de Kolmogorov-Smirnov¶

Test de Shapiro-Wilk¶

Proceso de Poisson¶

Ejercicio 5¶

Comparaciones con un GCL¶

Trabajo Práctico 2¶

Page Rank de Actores y Actrices¶

Ejercicio 2¶

Notas sobre el paper¶

El paper¶

El modelo¶

Ejercicio 3: Diseño del Web Service¶

Centralizado¶

Distribuido¶

Modificando los parametros¶

Media del tiempo de arribo¶

Lambda¶

Lambda1¶

Lambda2¶

p¶

Ejercicio 4¶

¿Recomienda a la entidad que implemente el cambio de cajero?¶

	jonvoight.html	angelinajolie.html	bradpitt.html	jenniferaniston.html	martinscorcese.html
jonvoight.html	0.0	1.0	1.0	0.0	0.0
angelinajolie.html	1.0	0.0	1.0	0.0	0.0
bradpitt.html	0.0	1.0	0.0	1.0	1.0
jenniferaniston.html	0.0	0.0	0.0	0.0	0.0
martinscorcese.html	0.0	0.0	0.0	0.0	0.0
robertdeniro.html	0.0	0.0	0.0	0.0	1.0